1 Contents/内容

以下のデータを整理

・The google mobility index. source:google mobility report(Accessed 25 July 2021)/https://www.google.com/covid19/mobility/

・The number of infected people and Coronavirus (COVID-19) Deaths/ COVID19 感染・死亡データ source:https://covid19.mhlw.go.jp/(Accessed 25 July 2021)

2 packages

rm(list = ls())

path <- getwd()

setwd(path)

# packages
pacman::p_load(tidyverse, plotly,readxl,DT, lubridate, extrafont)

# Font for windows and mac
if (stringr::str_detect(path, pattern="Users")){ 
  
   theme_set(theme_classic(base_size = 10, base_family = "HiraginoSans-W3"))  # For Mac OS
   base_family <- "HiraginoSans-W3"

 } else{
  
    theme_set(theme_classic(base_size = 10, base_family = "Arial")) # For Windows
   base_family <- "Arial"
}    

3 Read data

newly_confirmed_cases_daily <- readr::read_csv("input/newly_confirmed_cases_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))
confirmed_cases_cumulative_daily <-readr::read_csv("input/confirmed_cases_cumulative_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))

deaths_cumulative_daily <- readr::read_csv("input/deaths_cumulative_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))
JP_Mobility_Report <- readr::read_csv("input/2020_JP_Region_Mobility_Report.csv")
## Parsed with column specification:
## cols(
##   country_region_code = col_character(),
##   country_region = col_character(),
##   sub_region_1 = col_character(),
##   sub_region_2 = col_logical(),
##   metro_area = col_logical(),
##   iso_3166_2_code = col_character(),
##   census_fips_code = col_logical(),
##   place_id = col_character(),
##   date = col_date(format = ""),
##   retail_and_recreation_percent_change_from_baseline = col_double(),
##   grocery_and_pharmacy_percent_change_from_baseline = col_double(),
##   parks_percent_change_from_baseline = col_double(),
##   transit_stations_percent_change_from_baseline = col_double(),
##   workplaces_percent_change_from_baseline = col_double(),
##   residential_percent_change_from_baseline = col_double()
## )

4 load prefec_id

prefec_id <- read_excel("input/prefec_id.xlsx")

5 The number of cumulative confirmed cases/deaths at the end of June 2020/ COVID19 累積感染者数(6月30日)、累積死者(6月30日)

# Data on the number of cumulative deaths at the end of June 2020
cum_deaths_20200630 <- deaths_cumulative_daily %>% 
  mutate(year = year(Date),month=month(Date),day=day(Date)) %>% 
  filter(Date == "2020-06-30") %>% 
  mutate(.,month = 6,year=2020)

# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode

cum_deaths_20200630<- cum_deaths_20200630 %>% 
  mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))

# cumulative confirmed cases 
cum_confirmed_cases_20200630 <- confirmed_cases_cumulative_daily %>% 
  mutate(year = year(Date),month=month(Date),day=day(Date)) %>% 
  filter(Date == "2020-06-30") %>% 
  mutate(.,month = 6,year=2020)

# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode

cum_confirmed_cases_20200630<- cum_confirmed_cases_20200630 %>% 
  mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))

6 The number of cumulative confirmed cases/deaths at the end of May 2020/COVID19 累積感染者数(5月31日)、累積死者(5月31日)

# Data on the number of cumulative deaths at the end of may 2020
cum_deaths_20200531 <- deaths_cumulative_daily %>% 
  mutate(year = year(Date),month=month(Date),day=day(Date)) %>% 
  filter(Date == "2020-05-31") %>% 
  mutate(.,month = 5,year=2020)

# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode

cum_deaths_20200531<- cum_deaths_20200531 %>% 
  mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))



# cumulative confirmed cases 
cum_confirmed_cases_20200531 <- confirmed_cases_cumulative_daily %>% 
  mutate(year = year(Date),month=month(Date),day=day(Date)) %>% 
  filter(Date == "2020-05-31") %>% 
  mutate(.,month = 5,year=2020)

#  Recode "prefecture" to "prerec" and "All" to "zenkoku後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode

cum_confirmed_cases_20200531<- cum_confirmed_cases_20200531 %>% 
  mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))

7 moblity index for Japan

google_mobility_Japan2020  <- JP_Mobility_Report %>% 
  mutate(year = year(date),month=month(date),day=day(date)) %>% 
  # 2020年にデータを限定
  filter(year==2020) %>% 
  # 日本全国のsub_region1がnaなのでreplace_na
 replace_na(list(sub_region_1="zenkoku")) %>% 
  # sub_region1と月次で集計
  group_by(sub_region_1,month) %>% 
  summarise(
    retail_and_recreation = mean(retail_and_recreation_percent_change_from_baseline),
    grocery_and_pharmacy = mean(grocery_and_pharmacy_percent_change_from_baseline),
    parks = mean(parks_percent_change_from_baseline),
    transit_stations = mean(transit_stations_percent_change_from_baseline),
    workplaces = mean(workplaces_percent_change_from_baseline),
    residential = mean(residential_percent_change_from_baseline)
  )
## `summarise()` regrouping output by 'sub_region_1' (override with `.groups` argument)
  # make mobility index
  #  year=2020代入、sub_region1をprefecに変換。
google_mobility_Japan2020 <-google_mobility_Japan2020 %>% 
  mutate(.,prefec=sub_region_1,
         year=2020,
         google_mobility_index_4vari_average = (retail_and_recreation+grocery_and_pharmacy+transit_stations+workplaces)/4,
         mobility_index_6vari_average =(retail_and_recreation+grocery_and_pharmacy+transit_stations+workplaces+parks+residential)/6)

8 left_join 2020_may

covariates2020may <-left_join(cum_confirmed_cases_20200531, cum_deaths_20200531,by = c("prefec","year","month","day","Date","Prefecture"))

covariates2020may <- left_join(google_mobility_Japan2020, covariates2020may,by = c("prefec","month","year"))
covariates2020may <- left_join(covariates2020may,prefec_id,by="prefec")

9 left_join 2020_jun

covariates <-left_join(cum_confirmed_cases_20200630, cum_deaths_20200630,by = c("prefec","year","month","day","Date","Prefecture"))

covariates2020 <- left_join(google_mobility_Japan2020, covariates,by = c("prefec","month","year"))
covariates2020 <- left_join(covariates2020,prefec_id,by="prefec")

10 Save the data for May 2020./2020年5月を残して、csv保存

covariates202005 <- covariates2020may %>% 
# sub_region_1をドロップするために、ungroupを使う。googlemobliity indexの箇所でgroup_byを使っているため。
  filter(month==5) %>% 
  ungroup() %>% 
    select(-c(sub_region_1,Prefecture ))

covariates202005 <- covariates202005 %>% 
    rename(confirmed_cases_cumulative="Confirmed cases(Cumulative)",deaths_cumulative="Deaths(Cumulative)")

10.1 Write csv/データ保存

write_csv(covariates202005,"output/infection_death_mobility2020may.csv")

11 Save the data for June 2020./2020年6月を残して、csv保存

covariates202006 <- covariates2020 %>% 
  # sub_region_1をドロップするために、ungroupを使う。googlemobliity indexの箇所でgroup_byを使っているため。
  filter(month==6) %>% 
  ungroup() %>% 
    select(-c(sub_region_1,Prefecture ))
covariates202006 <- covariates202006 %>% 
    rename(confirmed_cases_cumulative="Confirmed cases(Cumulative)",deaths_cumulative="Deaths(Cumulative)")

11.1 write_csv/データ保存

write_csv(covariates202006,"output/infection_death_mobility2020Jun.csv")

12 Data check/データ確認

datatable(head(covariates2020))
summary(covariates2020)
##  sub_region_1           month    retail_and_recreation grocery_and_pharmacy
##  Length:528         Min.   : 2   Min.   :-52.419       Min.   :-10.767     
##  Class :character   1st Qu.: 4   1st Qu.:-11.797       1st Qu.: -0.500     
##  Mode  :character   Median : 7   Median : -5.525       Median :  1.403     
##                     Mean   : 7   Mean   : -7.880       Mean   :  1.346     
##                     3rd Qu.:10   3rd Qu.: -1.767       3rd Qu.:  3.363     
##                     Max.   :12   Max.   :  6.645       Max.   :  9.333     
##                                                                            
##      parks          transit_stations   workplaces        residential     
##  Min.   :-46.6129   Min.   :-58.16   Min.   :-42.3871   Min.   : 0.8387  
##  1st Qu.: -7.3473   1st Qu.:-27.40   1st Qu.:-14.8097   1st Qu.: 3.5667  
##  Median :  0.5333   Median :-19.84   Median :-10.3441   Median : 4.9000  
##  Mean   :  1.5555   Mean   :-20.45   Mean   :-11.1204   Mean   : 5.6391  
##  3rd Qu.:  8.8667   3rd Qu.:-11.25   3rd Qu.: -6.3169   3rd Qu.: 6.8065  
##  Max.   : 80.3226   Max.   : 15.57   Max.   : -0.3548   Max.   :21.4194  
##  NA's   :37         NA's   :2                                            
##     prefec               year      google_mobility_index_4vari_average
##  Length:528         Min.   :2020   Min.   :-39.048                    
##  Class :character   1st Qu.:2020   1st Qu.:-13.019                    
##  Mode  :character   Median :2020   Median : -8.440                    
##                     Mean   :2020   Mean   : -9.525                    
##                     3rd Qu.:2020   3rd Qu.: -4.445                    
##                     Max.   :2020   Max.   :  3.142                    
##                                    NA's   :2                          
##  mobility_index_6vari_average      Date             Prefecture       
##  Min.   :-25.995              Min.   :2020-06-30   Length:528        
##  1st Qu.: -8.739              1st Qu.:2020-06-30   Class :character  
##  Median : -5.317              Median :2020-06-30   Mode  :character  
##  Mean   : -5.330              Mean   :2020-06-30                     
##  3rd Qu.: -1.113              3rd Qu.:2020-06-30                     
##  Max.   : 14.952              Max.   :2020-06-30                     
##  NA's   :37                   NA's   :480                            
##  Confirmed cases(Cumulative)      day      Deaths(Cumulative)
##  Min.   :    0.00            Min.   :30    Min.   :  0.00    
##  1st Qu.:   43.75            1st Qu.:30    1st Qu.:  0.00    
##  Median :   82.00            Median :30    Median :  1.50    
##  Mean   :  763.31            Mean   :30    Mean   : 40.54    
##  3rd Qu.:  245.25            3rd Qu.:30    3rd Qu.: 19.75    
##  Max.   :18394.00            Max.   :30    Max.   :973.00    
##  NA's   :480                 NA's   :480   NA's   :480       
##  prefec_kanji             id       
##  Length:528         Min.   : 1.00  
##  Class :character   1st Qu.:12.75  
##  Mode  :character   Median :24.50  
##                     Mean   :25.56  
##                     3rd Qu.:36.25  
##                     Max.   :99.00  
## 
datatable(head(covariates202006))
summary(covariates202006)
##      month   retail_and_recreation grocery_and_pharmacy     parks        
##  Min.   :6   Min.   :-30.667       Min.   :-4.8667      Min.   :-33.767  
##  1st Qu.:6   1st Qu.:-10.442       1st Qu.:-0.1833      1st Qu.:-10.575  
##  Median :6   Median : -8.317       Median : 0.6667      Median : -3.517  
##  Mean   :6   Mean   : -9.025       Mean   : 0.9833      Mean   : -0.291  
##  3rd Qu.:6   3rd Qu.: -6.800       3rd Qu.: 2.1667      3rd Qu.:  9.100  
##  Max.   :6   Max.   :  2.867       Max.   : 5.4333      Max.   : 56.500  
##  transit_stations   workplaces       residential        prefec         
##  Min.   :-38.50   Min.   :-23.700   Min.   : 2.067   Length:48         
##  1st Qu.:-26.39   1st Qu.:-10.633   1st Qu.: 4.225   Class :character  
##  Median :-23.80   Median : -8.717   Median : 4.900   Mode  :character  
##  Mean   :-23.68   Mean   : -9.453   Mean   : 5.354                     
##  3rd Qu.:-19.57   3rd Qu.: -7.133   3rd Qu.: 5.983                     
##  Max.   :-12.40   Max.   : -5.433   Max.   :12.467                     
##       year      google_mobility_index_4vari_average
##  Min.   :2020   Min.   :-23.025                    
##  1st Qu.:2020   1st Qu.:-11.894                    
##  Median :2020   Median :-10.100                    
##  Mean   :2020   Mean   :-10.292                    
##  3rd Qu.:2020   3rd Qu.: -8.671                    
##  Max.   :2020   Max.   : -2.558                    
##  mobility_index_6vari_average      Date            confirmed_cases_cumulative
##  Min.   :-16.167              Min.   :2020-06-30   Min.   :    0.00          
##  1st Qu.: -8.375              1st Qu.:2020-06-30   1st Qu.:   43.75          
##  Median : -6.814              Median :2020-06-30   Median :   82.00          
##  Mean   : -6.018              Mean   :2020-06-30   Mean   :  763.31          
##  3rd Qu.: -4.300              3rd Qu.:2020-06-30   3rd Qu.:  245.25          
##  Max.   :  8.094              Max.   :2020-06-30   Max.   :18394.00          
##       day     deaths_cumulative prefec_kanji             id       
##  Min.   :30   Min.   :  0.00    Length:48          Min.   : 1.00  
##  1st Qu.:30   1st Qu.:  0.00    Class :character   1st Qu.:12.75  
##  Median :30   Median :  1.50    Mode  :character   Median :24.50  
##  Mean   :30   Mean   : 40.54                       Mean   :25.56  
##  3rd Qu.:30   3rd Qu.: 19.75                       3rd Qu.:36.25  
##  Max.   :30   Max.   :973.00                       Max.   :99.00

13 graph mobility index old

moblity_jp_zenkoku<-google_mobility_Japan2020 %>% 
  filter(., sub_region_1=="zenkoku") 


graph_mobility <- ggplot(data = moblity_jp_zenkoku, mapping = aes(x = month)) +
    theme_classic() +
      geom_line(aes(y= google_mobility_index_4vari_average, color = "Google mobility index(average of four mobility index)"))+
      geom_line(aes(y= mobility_index_6vari_average, color = "Google mobility index(average of six mobility index)"))+
        scale_x_continuous(breaks = seq(2,12,by=1)) +
   #凡例調整
      theme(legend.position="bottom")+
      theme(legend.background = element_rect(fill = NA))+
      theme(legend.direction  ="vertical")+
      labs(title="Google mobility index ",x = "month (in 2020)", y = "",color="") 

  ggplotly(graph_mobility)

14 save

if (stringr::str_detect(path, pattern="Users")){ 

  quartz(file = "output/graph_trends_graph_mobility_oldver.pdf", type = "pdf", 
         family = base_family, width = 7, height = 5,dpi=300, # for mac
        pointsize = 10)
print(graph_mobility)
dev.off()} else{

ggsave(graph_mobility,filename = "output/graph_trends_graph_mobility_oldver.pdf", 
       width = 7, height = 5,dpi=300) # for windows

}
## quartz_off_screen 
##                 2

15 graph mobility index

#2021Sep30 Waki

moblity_jp_zenkoku <- google_mobility_Japan2020 %>% 
  filter(., sub_region_1=="zenkoku") 

graph_mobility <- ggplot(data = moblity_jp_zenkoku, mapping = aes(x = month)) +
    theme_classic() +
      geom_line(aes(y= mobility_index_6vari_average, color = "Google mobility index(average of six mobility index)"))+
        scale_x_continuous(breaks = seq(2,12,by=1)) +
   #凡例調整
      theme(legend.position="none")+
      labs(title="",x = "Month in 2020", y = "",color="") +
      scale_color_manual(values = "black") +
      annotate("rect",xmin = 4,xmax = 5,ymin = -Inf,ymax = Inf, alpha = .3, fill = "gray") +
     annotate(geom = "text", x = 5, y = 0,
             label = "←The first state of emergency", size = 5, hjust = 0) + 
    theme(text = element_text(family = base_family))



ggplotly(graph_mobility)

16 save

#2021Sep30 Waki

if (stringr::str_detect(path, pattern="Users")){ 

  quartz(file = "output/graph_trends_graph_mobility.pdf", type = "pdf", 
         family = base_family, width = 7, height = 5,dpi=300, # for mac
        pointsize = 10)
print(graph_mobility)
dev.off()} else{

ggsave(graph_mobility,filename = "output/graph_trends_graph_mobility.pdf", 
       width = 7, height = 5,dpi=300) # for windows

}
## quartz_off_screen 
##                 2

17 Data checkデータ確認

datatable(moblity_jp_zenkoku)